library(tidyverse)
library(lubridate) # date manipulation
library(padr) # TS padding
library(zoo) # imputasi missing value TS
library(fpp) # TS dataset
library(TSstudio) # TS interactive viz
library(forecast) # algoritma forecasting
library(TTR) # SMA function
library(tseries) # adf.test
# supaya semua plot memiliki theme_minimal()
theme_set(theme_minimal())Time series data merupakan data yang terurut berdasarkan waktu dan disampel pada interval yang sama.
Forecasting merupakan suatu metode untuk memprediksi/meramalkan data di masa depan
Apa perbedaan Time Series Forecasting dengan Regresi?
Regresi adalah prediksi berdasarkan prediktor atau faktor yang mempengaruhinya. Time series forecasting memprediksi berdasarkan nilai data di masa lalu.
Regression
\[y = \beta_0+\beta_1*x_1+\beta_2*x_2+...+\beta_n*x_n\]
Time Series Forecasting
\[y_t = \beta_0+\beta_1*y_{t-1}+\beta_2*y_{t-2}+...+\beta_n*y_{t-n}\]
Ide utama dalam melakukan forecasting itu adalah korelasi dari data numerik, atau disebut sebagai autokorelasi.
Time series data: data yang berhubungan dengan waktu dan memiliki interval waktu yang tetap/sama.
Syarat utama data time series:
Knowledge Check
Apakah data berikut sudah memenuhi syarat data time series yang baik?
df_demand <- data.frame(
date = ymd(c("2021-5-3", "2021-5-4", "2021-5-6", "2021-5-7")),
demand = c(29, 79, 41, 88)
)
df_demandKesimpulan: Belum memenuhi karakteristik data TS yang baik, karena ada data yang terlongkap yaitu 5 Mei 2021.
df_price <- data.frame(
date = ymd(c("2021-5-16", "2021-5-19", "2021-5-18", "2021-5-17", "2021-5-20", "2021-5-21")),
price = c(1000, 1001, 1002, 1003, 1004, 1005)
)
df_priceKesimpulan: Belum terurut berdasarkan waktu (date)
Perbaikan yang dapat dilakukan sesuai syarat time series:
arrange()# mengurutkan df_price
df_price %>%
arrange(date)pad() dari package padrSecara default, pad() akan menambal tanggal berdasarkan kolom yang tipe datanya date. Mengisi nilai yang terlewat atau missing (NA), cara yang umum dilakukan dengan package zoo:
na.fill(): mengisi NA dengan sebuah nilai, Gunakan fill="extend" untuk mengisi dengan nilai rata-rata dengan nilai yang missingna.aggregate(): nilai aggregasi (mean, median)na.locf(): nilai terakhir sebelum missingNote: metode untuk mengisi missing value disesuaikan dengan perspektif dari businessnya.
# case: toko tutup di tanggal 5
# kita isi NA dengan nilai 0
df_demand %>%
pad() %>%
mutate(demand = na.fill(demand, 0))Kita akan menggunakan data emisi CO2 di Indonesia di mana datanya sudah tersimpan dalam folder data_input dengan nama environment_1970f.csv.
# read data
co2 <- read.csv("data_input/environment_1970f.csv")
head(co2)Data co2 terdiri dari 43 observasi yang mewakili kontribusi gas emisi per tahun terhadap atmosfer Indonesia (43 tahun, 1970-2012). Data ini terdiri dari 7 variabel, yaitu:
year: tahun.emisi co2 (kt): emisi yang berasal dari pembakaran bahan bakar fosil dan pembuatan semen, termasuk yang dihasilkan selama konsumsi.emisi co2 (metrik ton per kapita): idem.emisi metana (kt setara co2): emisi yang berasal dari aktivitas manusia (pertanian) dan dari produksi metana industri.emisi nitro oksida (ribu metrik ton setara co2): emisi dari pembakaran biomassa pertanian, kegiatan industri, dan pengelolaan ternak.emisi gas rumah kaca dan lainnya, HFC, PFC dan SF6 (ribu metrik ton setara co2): emisi hasil samping dari hidrofluorokarbon, perfluorokarbon, dan sulfur hexafluoride.total emisi gas rumah kaca (setara dengan co2): total CO2 tidak termasuk pembakaran biomassa siklus pendek (pembakaran limbah pertanian dan savannah), tetapi termasuk pembakaran biomassa lainnya (kebakaran hutan, pembusukan pasca-pembakaran, kebakaran gambut, dan pembusukan lahan gambut yang dikeringkan), semua sumber CH4 antropogenik, sumber N2O dan gas-F (HFC, PFC, dan SF6).Dari data co2 ini, kita akan menggunakan 2 kolom:
year untuk menunjukkan waktuCO2.emissions..metric.tons.per.capita. sebagai nilai yang kita amati untuk membuat object tsrange()range(co2$year)#> [1] 1970 2012
# memastikan apakah tahun sudah terurut dari 1970 sampai 2012
all(co2$year == seq(1970, 2012))#> [1] TRUE
seq(1970, 2012)#> [1] 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
#> [16] 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999
#> [31] 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012
Berikut adalah contoh cara penentuan frequency:
Knowledge Check
Tentukan frequency dari data time series berikut:
ts()Parameter:
data: vector numerik dari data time seriesstart: periode awal (opsional)
frequency: pola berulang pada data time seriesdata tahunan, kemudian kita mau lihat pola tahunan
# membuat object ts
co2_ts <- ts(data = co2$CO2.emissions..metric.tons.per.capita., # vektor numerik
start = 1970, # tahun awal
frequency = 1)
co2_ts#> Time Series:
#> Start = 1970
#> End = 2012
#> Frequency = 1
#> [1] 0.3119519 0.3306215 0.3580080 0.3954703 0.4021567 0.4128050 0.4612390
#> [8] 0.6002978 0.6677802 0.6601457 0.6426495 0.6634071 0.6822242 0.6640976
#> [15] 0.6944021 0.7347680 0.7229173 0.7184145 0.7552095 0.7348064 0.8243417
#> [22] 0.9735380 1.0788747 1.1452296 1.1416286 1.1420774 1.2669930 1.3738789
#> [29] 1.0218517 1.1599925 1.2452416 1.3748183 1.4102338 1.4364045 1.5098982
#> [36] 1.5084806 1.5015768 1.6118554 1.7638951 1.8651654 1.7679079 2.3335854
#> [43] 2.4089201
tsGunakan autoplot() dari package forecast, yang mengembalikan object ggplot:
library(forecast)
co2_ts %>%
autoplot()Misalkan ingin cek pola pada tahun 1990 - 2000 saja, gunakan window(start, end)
window: melakukan subsetting pada object time series
co2_ts %>%
window(start = 1990, end = 2000) %>%
autoplot()Insight: Pergerakan emisi gas co2 di Indonesia secara umum naik, tapi ada yang menarik di tahun 1997-1998. Terjadi peningkatan emisi gas yang amat tinggi pada tahun 1997 yang mungkin disebabkan oleh Kebakaran Hutan 1997
Diskusi: Objek co2_ts berasal dari data yang direkam tahunan dan kita mengatur pola tahunan (frequency = 1). Apakah kita bisa menganalisis pola berulang (musimannya)?
Jawaban: tidak bisa menganalisis pola berulang apabila frequency = 1
Note: pemilihan frequency umumnya menggunakan ukuran waktu 1 level di atasnya, atau yang lebih atas lagi. (data harian dengan seasonality mingguan / bulanan / dst.)
nybirth Datanybirth.csvbirth <- read.csv("data_input/nybirth.csv")
glimpse(birth)#> Rows: 168
#> Columns: 2
#> $ date <chr> "1946-01-01", "1946-02-01", "1946-03-01", "1946-04-01", "1946-0~
#> $ births <dbl> 26.663, 23.598, 26.931, 24.740, 25.806, 24.364, 24.477, 23.901,~
Data di atas merupakan data persentase kelahiran di New York per bulan. Terdiri dari:
date: tanggal saat dilakukan pencatatan persentase kelahiranbirths: persentase kelahiranCek class data birth:
class(birth)#> [1] "data.frame"
# masih data framedatedate untuk keperluan visualisasi nantinya# gunakan fungsi dari package lubridate
birth_clean <- birth %>%
mutate(date = ymd(date),
month = month(date, label = T))
birth_cleanrange(birth_clean$date)#> [1] "1946-01-01" "1959-12-01"
birth_clean %>%
arrange(date) %>%
pad() %>%
anyNA()#> [1] FALSE
Kesimpulan: data birth_clean sudah memenuhi ketiga karakteristik data TS, karena setelah diurutkan dan dilakukan padding, tidak ada missing value
Visualisasi data sebelum mengubah menjadi object ts untuk mengetahui pola berulang pada data birth
ggplot(data = birth_clean, aes(x = date, y = births)) +
geom_point() +
geom_line()ggplot(data = birth_clean, aes(x = date, y = births)) +
geom_point() +
geom_point(data = birth_clean %>%
filter(month %in% c("Jan", "Feb", "Jul")),
aes(col = month)) +
geom_line() +
scale_color_manual(values = c("red", "blue", "yellow"))Insight: Bagaimana pola dari data birth? ada pola yang berulang setiap 12 bulan (pola tahunan) karena setiap Feb tingkat kelahiran cenderung rendah, sedangkan tiap Jul cenderung tinggi.
birth_ts <- ts(data = birth_clean$births,
frequency = 12, # data bulanan, pola tahunan
start = c(1946, 1)) # opsional, data awal = 1946 Jan
birth_ts#> Jan Feb Mar Apr May Jun Jul Aug Sep Oct
#> 1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227
#> 1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110
#> 1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142
#> 1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907
#> 1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252
#> 1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110
#> 1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199
#> 1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462
#> 1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379
#> 1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784
#> 1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136
#> 1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484
#> 1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945
#> 1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012
#> Nov Dec
#> 1946 21.672 21.870
#> 1947 21.759 22.073
#> 1948 21.059 21.573
#> 1949 21.519 22.025
#> 1950 22.084 22.991
#> 1951 22.964 23.981
#> 1952 23.162 24.707
#> 1953 25.246 25.180
#> 1954 24.712 25.688
#> 1955 25.693 26.881
#> 1956 26.291 26.987
#> 1957 26.634 27.735
#> 1958 25.912 26.619
#> 1959 26.992 27.897
Misal kita tertarik untuk visualisasi birth dari Jan 1950 sampai Des 1955 dengan autoplot()
birth_ts %>%
window(start = c(1950, 1), end = c(1955, 12)) %>%
autoplot()Opsional: Gunakan ts_plot() dari package TSstudio untuk visualisasi interaktif:
birth_ts %>%
ts_plot()Decomposition adalah suatu tahapan dalam analisis time series untuk menguraikan data menjadi beberapa komponen dalam time series data, yaitu:
Untuk dapat menguraikan object time series kita menjadi 3 komponen tersebut, kita dapat menggunakan fungsi decompose().
birth_dc <- birth_ts %>%
decompose()Visualisasi hasil decompose:
birth_dc %>%
autoplot()Notes: Jika pada hasil decompose, trend masih membentuk sebuah pola maka dapat dicurigai masih ada seasonality yang belum ditangkap. Seharusnya trend cenderung naik atau cenderung turun secara smooth. Penyebabnya:
# coba buat ulang object ts dari birth_clean tapi frequency = 4
ts(data = birth_clean$births, frequency = 4) %>%
decompose() %>%
autoplot()Notes: Object time series dengan frequency = 1, tidak bisa dibuat decomposenya
# co2_ts %>%
# decompose()
# time series has no or less than 2 periodsTerdapat 2 jenis model pada data time series, yaitu:
\[Y_t = T_t + S_t + E_t\]
Data = Trend + Seasonal + Error
Contoh time series additive:
birth_ts %>%
autoplot()\[Y_t = T_t * S_t * E_t\]
Data = Trend * Seasonal * Error
Contoh time series multiplicative:
AirPassengers %>%
autoplot()Jika kita perhatikan lagi pada data birth_ts memiliki pola additive, karena varians dari polanya tetap atau konstan. Secara default, fungsi decompose() memiliki type = "additive".
Melakukan inspect komponen time series pada data birth_dc
1. Trend
Trend diperoleh dari hasil perhitungan center moving average (CMA). Tujuan utamanya untuk smoothing data sehingga diperoleh trend yang cenderung naik atau turun.
Berikut pendekatan manual untuk perhitungan trend, menggunakan fungsi ma(). Parameter:
x: object TSorder: berapa banyak data yang dilibatkan, menyesuaikan frequencycentre = T: menggunakan center moving average# perhitungan manual, di mana order = frequency pada object TS
ma(birth_ts, order = 12, centre = T) %>% head(24)#> Jan Feb Mar Apr May Jun Jul Aug
#> 1946 NA NA NA NA NA NA 23.98433 23.66213
#> 1947 22.35350 22.30871 22.30258 22.29479 22.29354 22.30562 22.33483 22.31167
#> Sep Oct Nov Dec
#> 1946 23.42333 23.16112 22.86425 22.54521
#> 1947 22.26279 22.25796 22.27767 22.35400
Bandingkan dengan $trend dari birth_dc:
# trend dari hasil decompose()
birth_dc$trend %>% head(24)#> Jan Feb Mar Apr May Jun Jul Aug
#> 1946 NA NA NA NA NA NA 23.98433 23.66213
#> 1947 22.35350 22.30871 22.30258 22.29479 22.29354 22.30562 22.33483 22.31167
#> Sep Oct Nov Dec
#> 1946 23.42333 23.16112 22.86425 22.54521
#> 1947 22.26279 22.25796 22.27767 22.35400
# hasilnya sama persis dengan perhitungan manual2. Seasonal
Additive: Data = Trend + Seasonal + Error Data de-trend: Data - Trend = Seasonal + Error
Data - Trend = Seasonal + Error
birth_detrend <- birth_ts - birth_dc$trendPendekatan manual
# mean tiap bulan
mean_month <- birth_detrend %>%
matrix(ncol = 12, byrow = T) %>%
colMeans(na.rm = T)
# mean global
mean_global <- mean(mean_month)
# nilai seasonality
mean_month - mean_global#> [1] -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
#> [7] 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
Bandingkan dengan $seasonal dari birth_dc:
# seasonal dari hasil decompose()
birth_dc$seasonal %>% head(12)#> Jan Feb Mar Apr May Jun
#> 1946 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
#> Jul Aug Sep Oct Nov Dec
#> 1946 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
# hasilnya sama persis dengan perhitungan manual3. Error
Additive: Data = Trend + Seasonal + Error Additive error: Error = Data - Trend - Seasonal
(birth_ts - birth_dc$trend - birth_dc$seasonal) %>%
autoplot()Bandingkan dengan $random dari birth_dc:
# error dari hasil decompose()
birth_dc$random %>%
autoplot()# hasilnya sama persis dengan perhitungan manualInti dari perhitungan manual:
Istilah:
$randomJika kita memiliki pola multiplicative pada data time series dan ingin membuat decomposenya cukup menambahkan parameter type = "multiplicative" pada fungsi decompose().
Ketika kita menemukan pola data kita mengandung multiplicative:
Cara 1: transformasi data multiplicative menjadi additive dengan fungsi
log. Setelah memperoleh hasil forecast kita dapat mengembalikan nilainya denganexp.
AirPassengers %>%
autoplot()# transformasi log
AirPassengers %>%
log() %>%
autoplot()(Opsional) Sifat logaritma: perkalian menjadi penjumlahan
\[y = T * S * E\] -> multiplicative \[log(y) = log(T * S * E)\] \[log(y) = log(T) + log(S) + log(E)\] -> additive
Cara 2: Tetap menggunakan model multiplicative, kemudian nanti hasil dibandingkan dengan memilih model dengan error yang paling kecil.
Parameter type dalam fungsi decompose(), secara default type = "additive"
air_decom <- AirPassengers %>%
decompose(type = "multiplicative")
air_decom %>%
autoplot()Trend
Menggunakan Center Moving Average (CMA)
ma(x = air_decom$x, order=12, centre = T) %>% head(12)#> Jan Feb Mar Apr May Jun Jul Aug
#> 1949 NA NA NA NA NA NA 126.7917 127.2500
#> Sep Oct Nov Dec
#> 1949 127.9583 128.5833 129.0000 129.7500
Bandingkan dengan trend di air_decom:
air_decom$trend %>% head(12)#> Jan Feb Mar Apr May Jun Jul Aug
#> 1949 NA NA NA NA NA NA 126.7917 127.2500
#> Sep Oct Nov Dec
#> 1949 127.9583 128.5833 129.0000 129.7500
# hasilnya sama persis dengan perhitungan manualSeasonality
Multiplicative: Data = Trend * Seasonal * Error De-trend: Data / Trend = Seasonal * Error
seasxerr <- air_decom$x /air_decom$trend
# mean of each month
mean_month <- seasxerr %>%
matrix(ncol = 12, byrow = T) %>%
colMeans(na.rm = T)
# mean global
mean_global <- mean(mean_month)
# Seasonality Calculation
mean_month / mean_global#> [1] 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
#> [8] 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
Bandingkan dengan seasonal di air_decom:
air_decom$seasonal %>% head(12)#> Jan Feb Mar Apr May Jun Jul
#> 1949 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
#> Aug Sep Oct Nov Dec
#> 1949 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244
# hasilnya sama persis dengan perhitungan manualError
Multiplicative: Data = Trend * Seasonal * Error Multiplicative Error: Error = Data / Trend / Seasonal
(air_decom$x/air_decom$trend/air_decom$seasonal) %>%
autoplot()Bandingkan dengan random di air_decom:
air_decom$random %>%
autoplot()# hasilnya sama persis dengan perhitungan manualSeasonality analysis membantu kita mengetahui di waktu mana saja yang nilai datanya tinggi/rendah pada periode seasonal yang kita amati. Kita bisa menggunakan informasi $seasonal untuk membuat visualisasi seasonal bar plot:
# analisis seasonal dengan bar plot
birth_clean %>%
mutate(seasonal = birth_dc$seasonal) %>% # mengekstrak nilai seasonal
distinct(month, seasonal) %>% # mengambil baris yang unique saja
ggplot(aes(x = month, y = seasonal)) + # melakukan visualisasi
geom_col()Seasonality adjusted adalah data time series yang sudah dibuang efek seasonal nya. Umumnya digunakan untuk lebih mudah mendeteksi error/kejadian luar biasa/anomali dari data (tidak terganggu efek seasonal). Hal ini untuk kebutuhan exploratory data saja.
Berikut contoh data birth_ts yang sudah dibuang efek seasonalnya:
Model additive:
Data = Trend + Seasonality + Error Data - Seasonality = Trend + Error
# $x = data asli
# $seasonal = seasonality
(birth_dc$x - birth_dc$seasonal) %>%
autoplot() %>%
plotly::ggplotly()Visualisasi berikut menunjukkan line plot untuk Trend sekaligus Seasonal Adjusted data. Dengan seperti ini, kita dapat dengan mudah mendeteksi anomali (titik-titik data yang jauh dari trend nya).
Credit to Pak Bayu Andrianto :)
birth_dc$trend %>%
autoplot(series = "Trend") +
autolayer(birth_dc$x - birth_dc$seasonal, series = "Seasonal adjusted")ts(). Parameter:data: vektor numerikstart: tanggal awalfrequency: pola berulang yang ingin ditangkapautoplot():decompose() menjadi 3 komponen:Tujuan melakukan decomposition: - insight dari trend dan seasonal - mengetahui apakah frequency yang kita set sudah tepat atau belum.
DAY 2
Metode:
Metode ini menggunakan rataan bergerak (moving average) untuk melakukan smoothing data. Karena menggunakan rataan, maka bobot yang digunakan sama untuk setiap observasi di masa lalu.
\[y_t = \frac{y_{t-1}+y_{t-2}+...+y_{t-n}}{n}\]
SMA tepat digunakan untuk data yang tidak mengandung trend dan tidak ada seasonal
trend: pola data secara umum apakah naik/turun seasonal: pola berulang
y <- c(3, 4, 5, 6)
mean(y)#> [1] 4.5
(3+4+5+6)/4#> [1] 4.5
1/4 * 3 + 1/4 * 4 + 1/4 * 5 + 1/4 * 6#> [1] 4.5
SMA: berdasarkan rata-rata, dan tiap datanya memiliki bobot yang sama
Kita menggunakan data curah hujan tahunan dari tahun 1813-1912 (100 tahun) di London, England.
library(TTR)
rain <- scan("data_input/precip1.dat", skip = 1)
head(rain)#> [1] 23.56 26.07 21.86 31.24 23.65 23.88
Note: skip = 1 untuk melongkap 1 observasi pada file .dat, karena ada tulisan “Total”.
data tahunan -> pola tahunan
rain_ts <- ts(data = rain,
start = 1813, # tahun awal
frequency = 1)
rain_ts#> Time Series:
#> Start = 1813
#> End = 1912
#> Frequency = 1
#> [1] 23.56 26.07 21.86 31.24 23.65 23.88 26.41 22.67 31.69 23.86 24.11 32.43
#> [13] 23.26 22.57 23.00 27.88 25.32 25.08 27.76 19.82 24.78 20.12 24.34 27.42
#> [25] 19.44 21.63 27.49 19.43 31.13 23.09 25.85 22.65 22.75 26.36 17.70 29.81
#> [37] 22.93 19.22 20.63 35.34 25.89 18.65 23.06 22.21 22.18 18.77 28.21 32.24
#> [49] 22.27 27.57 21.59 16.93 29.48 31.60 26.25 23.40 25.42 21.32 25.02 33.86
#> [61] 22.67 18.82 28.44 26.16 28.17 34.08 33.82 30.28 27.92 27.14 24.40 20.35
#> [73] 26.64 27.01 19.21 27.74 23.85 21.23 28.15 22.61 19.80 27.94 21.47 23.52
#> [85] 22.86 17.69 22.54 23.28 22.17 20.84 38.10 20.65 22.97 24.26 23.01 23.67
#> [97] 26.75 25.36 24.79 27.88
Apakah data rain dapat di-decompose? Karena frequency = 1, tidak bisa dilakukan decompose
rain_ts %>%
autoplot()Istilah: Data curah hujan tidak memiliki trend dan seasonal, atau bisa kita sebut sebagai data stasioner.
Fungsi yang digunakan untuk fitting model dengan metode SMA adalah SMA() dari package TTR dengan parameter:
x: objek time seriesn: ordorain_sma3 <- SMA(rain_ts, n = 3) # melihat 3 data ke belakang
rain_sma3#> Time Series:
#> Start = 1813
#> End = 1912
#> Frequency = 1
#> [1] NA NA 23.83000 26.39000 25.58333 26.25667 24.64667 24.32000
#> [9] 26.92333 26.07333 26.55333 26.80000 26.60000 26.08667 22.94333 24.48333
#> [17] 25.40000 26.09333 26.05333 24.22000 24.12000 21.57333 23.08000 23.96000
#> [25] 23.73333 22.83000 22.85333 22.85000 26.01667 24.55000 26.69000 23.86333
#> [33] 23.75000 23.92000 22.27000 24.62333 23.48000 23.98667 20.92667 25.06333
#> [41] 27.28667 26.62667 22.53333 21.30667 22.48333 21.05333 23.05333 26.40667
#> [49] 27.57333 27.36000 23.81000 22.03000 22.66667 26.00333 29.11000 27.08333
#> [57] 25.02333 23.38000 23.92000 26.73333 27.18333 25.11667 23.31000 24.47333
#> [65] 27.59000 29.47000 32.02333 32.72667 30.67333 28.44667 26.48667 23.96333
#> [73] 23.79667 24.66667 24.28667 24.65333 23.60000 24.27333 24.41000 23.99667
#> [81] 23.52000 23.45000 23.07000 24.31000 22.61667 21.35667 21.03000 21.17000
#> [89] 22.66333 22.09667 27.03667 26.53000 27.24000 22.62667 23.41333 23.64667
#> [97] 24.47667 25.26000 25.63333 26.01000
autolayer()rain_ts %>%
autoplot(series = "Rain") +
autolayer(rain_sma3, series = "SMA(3)")# series = label di legendKnowledge Check: Misal kita gunakan SMA3 untuk forecast tahun 1913, maka nilainya adalah 26.01
# manual
tail(rain_ts, 3) %>% mean()#> [1] 26.01
# prediksi nilai di masa depan
library(forecast)
forecast(rain_sma3, h = 1)#> Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
#> 1913 26.00996 23.68437 28.33555 22.45328 29.56665
Kelemahan SMA:
n data ke belakangPada exponential smoothing, bobot yang digunakan menurun secara eksponensial. Data yang lebih baru akan diberikan bobot yang lebih besar dibandingkan data yang sudah lama. Metode exponential smoothing terbagi menjadi 3 yaitu:
Fungsi yang digunakan untuk membuat model exponential smoothing, yaitu:
HoltWinters(), parameter yang digunakan, yaitu:x: objek time series
secara default parameter alpha, beta, dan gamma adalah NULL. Parameter harus diubah menjadi FALSE apabila komponen tersebut tidak ada.
alpha = smoothing untuk error, nilainya antara 0 - 1beta = smoothing untuk trend, nilainya antara 0 - 1gamma = smoothing untuk seasonal, nilainya antara 0 - 1seasonality = jenis model, apakah additive/multiplicativeets -> bobot abg
SES: HoltWinters(x, beta = F, gamma = F)
\[y_{t+1} = \alpha y_{t} + \alpha (1-\alpha)y_{t-1}+\alpha(1-\alpha)^2y_{t-2}+...\]
Double/Holt Exponential Smoothing: HoltWinters(x, gamma = F)
\[y_{t+h|t} = \ell_t +hb_t\] \[\ell_t = \alpha y_t+(1-\alpha)(\ell_{t-1}+b_{t-1})\] \[b_t = \beta^{*}(\ell-\ell_{t-1}) +(1-\beta^{*})b_{t-1}\]
Triple/Holt Winters Exponential Smoothing: HoltWinters(x)
\[y_{t+h|t} = \ell_t +hb_t+s_{t+h-m(k+1)}\] \[\ell_t = \alpha y_t+(1-\alpha)(\ell_{t-1}+b_{t-1})\] \[b_t = \beta^{*}(\ell-\ell_{t-1}) +(1-\beta^{*})b_{t-1}\] \[s_t = \gamma(y_t-\ell_{t-1}-b_{t-1})+(1-\gamma)s_{t-m}\]
SES baik digunakan untuk data time series yang tidak mengandung trend dan tidak seasonal. SES akan melakukan smoothing terhadap komponen error/random.
Kita akan menggunakan data curah hujan untuk melakukan fitting model dan forecasting menggunakan metode SES.
rain_ts %>%
autoplot()Cross validation: train-test splitting
# 20 tahun terakhir untuk test
rain_test <- tail(rain_ts, 20)
# 80 tahun untuk train (sisanya)
rain_train <- head(rain_ts, -20) # ekuivalen dengan head(rain_ts, 80)Fitting Model dengan fungsi HoltWinters()
ets => parameter abg
no trend, no seasonal -> beta = F, gamma = F
rain_ses <- HoltWinters(rain_train, beta = F, gamma = F)
rain_ses#> Holt-Winters exponential smoothing without trend and without seasonal component.
#>
#> Call:
#> HoltWinters(x = rain_train, beta = F, gamma = F)
#>
#> Smoothing parameters:
#> alpha: 0.03595009
#> beta : FALSE
#> gamma: FALSE
#>
#> Coefficients:
#> [,1]
#> a 25.24619
Note: secara default fungsi HoltWinters() akan mencari parameter smoothing yang menurutnya optimal -> meminimalisir error di data train. “Optimal” itu belum tentu menghasilkan the best performance untuk di data test.
Apabila ingin set nilai smoothing sesuai keinginan, misal alpha di nilai 0.8 maka: HoltWinters(rain_train, alpha = 0.8, beta = F, gamma = F)
Forecasting dengan fungsi forecast()
rain_forecast <- forecast(rain_ses, h = 20) # panjang data test
rain_forecast %>%
autoplot()Parameter: h = horizon (berapa waktu kedepan yang ingin di-forecast?)
Output:
$mean): nilai forecasting di masa depanVisualisasi data actual dengan hasil forecasting SES
rain_ts %>%
autoplot(series = "Actual") +
autolayer(rain_forecast$fitted, series = "Train") +
autolayer(rain_forecast$mean, series = "Test")Note:
$fitted untuk hasil forecast di data train$mean untuk hasil forecast di data testEvaluasi model dengan fungsi accuracy() dari library forecast
Parameter: object_forecast dan object_ts_test
accuracy(rain_forecast, rain_test)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.5937162 4.310376 3.431055 -0.4116781 13.7016 0.6794168
#> Test set -1.2686879 4.280286 3.248594 -7.8305099 13.7201 0.6432859
#> ACF1 Theil's U
#> Training set -0.05786692 NA
#> Test set -0.24464714 0.7232343
Train: 13.7% Test: 13.72%
Metric evaluasi:
Holt digunakan untuk data time series yang mengandung trend, namun tidak seasonal. Holt akan melakukan smoothing terhadap komponen error/random dan trend.
Kita akan menggunakan data kontribusi emisi gas CO2 di Indonesia untuk melakukan fitting model dan forecasting menggunakan metode holt.
co2_ts %>%
autoplot()Note: ingat, karena frequency = 1 kita tidak dapat melakukan decompose()
Cross validation
# 3 tahun terakhir untuk test
co2_test <- tail(co2_ts, 3)
# 40 tahun untuk train (sisanya)
co2_train <- head(co2_ts, -3)Fitting Model Double Exponential Smoothing (DES)
ets => parameter abg
ada trend, tidak ada seasonal -> gamma = F
co2_des <- HoltWinters(co2_train, gamma = F)
co2_des#> Holt-Winters exponential smoothing with trend and without seasonal component.
#>
#> Call:
#> HoltWinters(x = co2_train, gamma = F)
#>
#> Smoothing parameters:
#> alpha: 0.9035708
#> beta : 0.03548442
#> gamma: FALSE
#>
#> Coefficients:
#> [,1]
#> a 1.85771545
#> b 0.03875735
Hati-hati: Jangan set nilai parameter smoothing menjadi TRUE. Contoh kalau beta = T artinya beta = 1
Forecasting
co2_forecast <- forecast(co2_des, h = 3) # panjang data testVisualisasi data actual dengan hasil forecasting
co2_ts %>%
autoplot(series = "Actual") +
autolayer(co2_forecast$fitted, series = "Train") +
autolayer(co2_forecast$mean, series = "Test")Evaluasi model
forecast::accuracy(co2_forecast, co2_test)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.01648723 0.08542822 0.05645794 1.323059 5.813989 0.9019799
#> Test set 0.23490765 0.34851219 0.32061760 9.284487 14.132587 5.1222313
#> ACF1 Theil's U
#> Training set -0.01284696 NA
#> Test set -0.13435655 0.9092797
MAPE:
Tindak lanjut:
ts_grid() pada package TSstudioHolt winters digunakan untuk data time series yang mengandung trend dan seasonal. Holt winters melakukan smoothing terhadap komponen error/random (alpha), trend (beta), dan seasonal (gamma). Kita akan menggunakan data persentase kelahiran di New York untuk melakukan fitting model dan forecasting menggunakan metode holt winters.
birth_ts %>%
autoplot()Cross validation
# 3 tahun (36 bulan) terakhir untuk test
birth_test <- tail(birth_ts, 3*12)
# sisanya untuk train
birth_train <- head(birth_ts, -3*12)Fitting Model
Secara default HoltWinters akan menganggap data sebagai additive seasonality (seasonal = "additive"). Apabila data berupa multiplicative, gunakan parameter seasonal = "multiplicative" di dalam fungsi Holtwinters().
birth_hw <- HoltWinters(birth_train, seasonal = "additive")
birth_hw#> Holt-Winters exponential smoothing with trend and additive seasonal component.
#>
#> Call:
#> HoltWinters(x = birth_train, seasonal = "additive")
#>
#> Smoothing parameters:
#> alpha: 0.8238483
#> beta : 0.01879443
#> gamma: 1
#>
#> Coefficients:
#> [,1]
#> a 27.36749822
#> b 0.02551976
#> s1 -0.62729995
#> s2 -1.71011964
#> s3 1.64280798
#> s4 0.05574370
#> s5 1.08739310
#> s6 0.27050131
#> s7 1.54289802
#> s8 0.92194413
#> s9 0.29164988
#> s10 0.48480479
#> s11 -1.50609008
#> s12 -0.38049822
Forecasting
birth_forecast <- forecast(birth_hw, h = 3*12)horizon = kita mau forecast berapa data ke depan?
Visualisasikan data actual dengan hasil forecasting holt winters
birth_ts %>%
autoplot(series = "Actual") +
autolayer(birth_forecast$fitted, series = "Train") +
autolayer(birth_forecast$mean, series = "Test")Evaluasi Model
forecast::accuracy(birth_forecast, birth_test)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.09810363 0.7530131 0.5681970 0.3510623 2.347179 0.5687989
#> Test set -0.52267504 1.0812812 0.8125558 -2.0048640 3.015850 0.8134166
#> ACF1 Theil's U
#> Training set 0.1490922 NA
#> Test set 0.6856166 0.6939018
MAPE:
ets() dari library forecast yang merupakan error, trend, dan seasonal. Parameter yang digunakan, yaitu:y: objek time series.model: model/pola time series untuk error, trend, dan seasonal.
A: additiveM: multiplicativeN: noneZ: auto -> berdasarkan AIC (information loss, semakin kecil semakin baik)alpha: bobot untuk smoothing error, nilainya antara 0 - 1.beta: bobot untuk smoothing trend, nilainya antara 0 - 1.gamma: bobot untuk smoothing seasonal, nilainya antara 0 - 1.Jika bobot semakin mendekati 1, maka nilai prediksi/forecast lebih mempertimbangkan data terbaru.
Single Exponential Smooting/SES: ets(y, model = "*NN")
Double Exponential Smooting/Holt: ets(y, model = "**N")
Triple Exponential Smooting/Holt Winters: ets(y, model = "***")
birth_ts %>%
autoplot()Trend: cenderung linear -> additive Seasonality: varians konstan -> additive
birth_ets <- ets(birth_train, model = "AAA")
birth_ets#> ETS(A,Ad,A)
#>
#> Call:
#> ets(y = birth_train, model = "AAA")
#>
#> Smoothing parameters:
#> alpha = 0.7193
#> beta = 0.0001
#> gamma = 0.0001
#> phi = 0.863
#>
#> Initial states:
#> l = 26.3085
#> b = -0.6758
#> s = -0.5129 -1.2048 0.7337 0.6324 1.1251 1.526
#> -0.0403 0.3452 -0.7701 0.9206 -2.1293 -0.6257
#>
#> sigma: 0.6392
#>
#> AIC AICc BIC
#> 544.1604 550.2135 596.0508
Cobalah dengan parameter model:
Forecasting
birth_ets_forecast <- forecast(birth_ets, h = 3*12) # panjang data testVisualisasikan data actual dengan hasil forecasting holt winters
birth_ts %>%
autoplot(series = "Actual") +
autolayer(birth_ets_forecast$fitted, series = "Train") +
autolayer(birth_ets_forecast$mean, series = "Test")Evaluasi Model
forecast::accuracy(birth_ets_forecast, birth_test)#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.05781606 0.5965746 0.4476847 0.1980528 1.827651 0.4481590
#> Test set -0.05731999 0.7280613 0.5871715 -0.2885430 2.142523 0.5877936
#> ACF1 Theil's U
#> Training set 0.1066730 NA
#> Test set 0.5234464 0.4572179
Bandingkan performa model forecasting data birth dengan fungsi HoltWinters() dengan ets():
Fungsi HoltWinters:
Fungsi ets (AAA):
DAY 3
Auto Regressive Integrated Moving Average (ARIMA) adalah gabungan dari dua metode, yaitu Autoregressive (AR) dan Moving Average (MA). I-nya menjelaskan Integrated. Konsep utama dari ARIMA, yaitu menjelaskan autokorelasi pada data.
AR(p): Autoregressive
Pada model regresi linier, target (\(y\)) diprediksi berdasarkan prediktor (\(x\)).
\[\hat{y}=\beta_0+\beta_1*x_1+...+\beta_n*x_n\]
Pada model autoregressive, target di masa depan (\(\hat{y_t}\)) diforecast berdasarkan nilai target di masa lalu (\(y_{t-1}, ..., y_{t-n}\)).
\[\hat{y}_t = c + \beta_1.y_{t-1} + \cdots + \beta_p.y_{t-p} + e_{t}\]
Dikenal istilah lag/delay, yaitu data pada waktu sebelumnya, di mana:
Contoh pada data harian:
Order \(p\): berapa banyak lag untuk nilai \(y\) yang kita gunakan dalam melakukan forecasting
MA(q): Moving average
Nilai \(\hat{y_t}\) dihasilkan dari smoothing error/random. Hasil smoothing yang dilakukan pada error/random bukan digunakan untuk melakukan forecast nilai \(\hat{y_t}\), melainkan digunakan untuk mengestimasi trend.
\[\hat{y}_t = c + \theta_1.\varepsilon_{t-1} + \cdots + \theta_q.\varepsilon_{t-q} + e_{t}\]
Order \(q\): berapa banyak lag untuk nilai error \(\varepsilon\) yang kita gunakan dalam melakukan forecasting
ARIMA(p,d,q): non-seasonal ARIMA
\[\hat{y}_t = c + \beta_1.y_{t-1} + \cdots + \beta_p.y_{t-p} + \theta_1.\varepsilon_{t-1} + \cdots + \theta_q.\varepsilon_{t-q} + e_{t}\]
Syarat ARIMA: time series harus stasioner (rata-rata dan variansinya konstan), artinya tidak ada trend dan tidak ada seasonal. Jika data time series belum stationer, maka perlu dilakukan differencing/Integrated pada ARIMA.
Cara yang paling umum untuk membuat data time series menjadi stasioner adalah dengan melakukan differencing diff(), yaitu mengurangi data saat ini dengan data sebelumnya (lag = 1). Terkadang, tergantung pada kompleksitas data, jumlah differencing bisa lebih dari 1 kali.
x <- c(5, 9, 3, 4, 7)
x %>%
diff()#> [1] 4 -6 1 3
Pertanyaannya, mengapa kita perlu melakukan differencing untuk mendapatkan time series yang stasioner? Berikut adalah penjelasannya:
Data AirPassengers sebelum differencing memiliki korelasi sangat kuat antar lag nya:
data.frame(
xt = c(AirPassengers),
lag1 = lag(x = as.vector(AirPassengers), n = 1),
lag2 = lag(x = as.vector(AirPassengers), n = 2),
lag3 = lag(x = as.vector(AirPassengers), n = 3),
lag4 = lag(x = as.vector(AirPassengers), n = 4),
lag5 = lag(x = as.vector(AirPassengers), n = 5)
) %>%
na.omit() %>%
GGally::ggcorr(label = T)Data AirPassengers setelah differencing, korelasi antar lag tidak terlalu kuat.
data.frame(
xt = c(diff(AirPassengers)),
lag1 = lag(x = as.vector(diff(AirPassengers)), n = 1),
lag2 = lag(x = as.vector(diff(AirPassengers)), n = 2),
lag3 = lag(x = as.vector(diff(AirPassengers)), n = 3),
lag4 = lag(x = as.vector(diff(AirPassengers)), n = 4),
lag5 = lag(x = as.vector(diff(AirPassengers)), n = 5)
) %>%
na.omit() %>%
GGally::ggcorr(label = T)Ide utama melakukan differencing adalah ketika melakukan prediksi, agar multicolinearity antar lag-nya tidak terlalu besar.
no <- read.csv("data_input/environment_1970f.csv")Data terdiri dari 43 observasi yang mewakili kontribusi gas emisi per tahun terhadap atmosfer Indonesia (43 tahun, 1970-2012). Sebelumnya kita telah menggunakan kolom emisi CO2, sekarang mari kita gunakan informasi emisi gas Nitro Oksida (ribu metrik ton setara CO2) yaitu emisi dari pembakaran biomassa pertanian, kegiatan industri, dan pengelolaan ternak.
p <- no %>%
ggplot(aes(x = year,
y = Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent.)) +
geom_point() +
geom_line() +
labs(y = "Nitrous Oxide")
plotly::ggplotly(p) # interactive plottingBerdasarkan line plot, apakah dapat dikatakan data sudah stasioner?
Dari visualisasi kita belum yakin secara objektif apakah data kita memiliki trend/seasonal.
no_ts <- ts(data = no$Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent.,
start = 1970,
frequency = 1)# 8 tahun terakhir untuk data test, sisanya jadikan data train
no_test <- tail(no_ts, 8)
no_train <- head(no_ts, -8)ADF test
Kalau hanya dari visualisasi, penentuan apakah time series sudah stasioner akan sangat subjektif. Untuk itu, kita dapat melakukan Augmented Dickey-Fuller (ADF) test untuk menguji stationarity dengan fungsi adf.test() dari library tseries.
Hipotesis:
Harapannya adalah p-value < alpha (0.05) -> stasioner
Note: lakukan tes hanya pada data train saja
library(tseries)
adf.test(no_train)#>
#> Augmented Dickey-Fuller Test
#>
#> data: no_train
#> Dickey-Fuller = -3.2514, Lag order = 3, p-value = 0.09523
#> alternative hypothesis: stationary
Berdasarkan ADF test, data no_train sudah stasioner atau belum?
p-value = 0.09523 alpha = 0.05
p-value > alpha, kesimpulannya gagal tolak H0 sehingga tidak stasioner
Apabila data belum stasioner, maka lakukan differencing hingga stasioner.
no_train_diff <- no_train %>% diff()Uji lagi ADF test: harapan p-value < alpha (0.05)
no_train_diff %>%
adf.test()#>
#> Augmented Dickey-Fuller Test
#>
#> data: .
#> Dickey-Fuller = -5.943, Lag order = 3, p-value = 0.01
#> alternative hypothesis: stationary
p-value = 0.01 < alpha (0.05) sehingga tolak H0
Kesimpulan: Data no_train perlu dilakukan differencing sebanyak satu kali hingga stasioner.
Mari kita bandingkan no_train tanpa dan dengan differencing melalui visualisasi:
no_train %>%
autoplot(series = "Actual") +
autolayer(no_train_diff, series = "Diff 1")ARIMA(p,d,q) = AR(p) I(d) MA(q)
Banyaknya differencing yang dilakukan sebelumnya menjadi bagian order d untuk komponen Integrated I(d), sehingga menjadi model ARIMA(p,1,q)
auto.arima()PENTING: data yang di-fit pada fungsi auto.arima() adalah data sebelum differencing
Gunakan parameter seasonal = F untuk ARIMA biasa
no_auto <- auto.arima(y = no_train, seasonal = F)
summary(no_auto)#> Series: no_train
#> ARIMA(0,1,1)
#>
#> Coefficients:
#> ma1
#> -0.8150
#> s.e. 0.0919
#>
#> sigma^2 estimated as 3317659997: log likelihood=-420.96
#> AIC=845.93 AICc=846.32 BIC=848.98
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE ACF1
#> Training set 10972.01 55929.24 31826.05 0.950838 23.63757 0.7418046 -0.04361855
Apa artinya ARIMA(0,1,1)?
Arima() dari package forecast# membuat model dengan order p dan q nya asal
# misal ARIMA(1,1,2)
no_manual <- Arima(y = no_train, order = c(1,1,2))
summary(no_manual)#> Series: no_train
#> ARIMA(1,1,2)
#>
#> Coefficients:
#> ar1 ma1 ma2
#> -0.8050 0.1723 -0.8276
#> s.e. 0.1183 0.1683 0.1480
#>
#> sigma^2 estimated as 3169988254: log likelihood=-419.85
#> AIC=847.7 AICc=849.08 BIC=853.8
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE ACF1
#> Training set 10336.83 52987.77 29322.74 1.500105 21.51614 0.6834571 -0.119463
Bandingkan kedua performa model:
Silahkan jadikan model auto sebuah baseline model, tapi perlu ingat bahwa model auto belum tentu menjadi model yang terbaik
Menentukan ARIMA (p, d, q) secara manual
Langkah 1: Tentukan d (dari differencing)
Langkah 2: Tentukan p dan q
tribble(~model, ~ACF, ~PACF,
"AR(p)", "Die Down", "Cut off after lag p",
"MA(q)", "Cut off after lag q", "Die Down",
"ARMA(p, q)", "Die Down", "Die Down",
"No AR(p) or MA(q)", "No spike", "No spike"
)Perbedaan PACF dan ACF:
Misal hari ini Kamis:
Die down vs Cut off
Mari kita cek ACF dan PACF untuk no_train_diff
PENTING: Saat plot ACF dan PACF, gunakan data yang sudah di-differencing (karena kita akan melakukannya secara manual)
# fungsi tsdisplay dari package forecast
tsdisplay(no_train_diff)Melakukan plot ACF dan PACF bisa menggunakan tsdisplay() namun agar tampil lebih besar, kita buat custom function plot_cf() sebagai berikut:
# opsional: custom function
library(patchwork)
plot_cf <- function(x, lag.max = NULL){
p1 <- ggAcf(x, lag.max=lag.max) + labs(title = "Plot ACF and PACF")
p2 <- ggPacf(x, lag.max=lag.max) + labs(title = NULL)
p1 / p2
}plot_cf(no_train_diff)Garis putus-putus horizontal adalah batasan signifikansi autokorelasi, range nilainya adalah +- 1.96 / akar(T), dengan T adalah panjang time series.
Dari plot ACF-PACF dapat disimpulkan:
Setelah itu mari tentukan order:
Kombinasi model ARIMA(p,d,q):
Istilah: Spike -> lag dengan nilai autokorelasi yang signifikan
Fitting kombinasi model
Note: fit data yang sebelum differencing
no_arima1 <- Arima(y = no_train, order = c(1,1,1))
no_arima2 <- Arima(y = no_train, order = c(2,1,1))
no_arima3 <- Arima(y = no_train, order = c(4,1,1))Cek error di training set dengan summary()
summary(no_arima1)#> Series: no_train
#> ARIMA(1,1,1)
#>
#> Coefficients:
#> ar1 ma1
#> 0.0267 -0.8216
#> s.e. 0.1943 0.1019
#>
#> sigma^2 estimated as 3420512259: log likelihood=-420.96
#> AIC=847.91 AICc=848.71 BIC=852.49
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE ACF1
#> Training set 11073.59 55922.5 31760.79 1.07617 23.49991 0.7402834 -0.06073625
summary(no_arima2)#> Series: no_train
#> ARIMA(2,1,1)
#>
#> Coefficients:
#> ar1 ar2 ma1
#> -0.0445 -0.2210 -0.7471
#> s.e. 0.2000 0.1822 0.1416
#>
#> sigma^2 estimated as 3379874911: log likelihood=-420.28
#> AIC=848.55 AICc=849.93 BIC=854.66
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE
#> Training set 10008.72 54713.83 30326.12 0.6124584 22.41164 0.7068441
#> ACF1
#> Training set -0.06979143
summary(no_arima3)#> Series: no_train
#> ARIMA(4,1,1)
#>
#> Coefficients:
#> ar1 ar2 ar3 ar4 ma1
#> -0.5841 -0.6503 -0.4078 -0.3948 -0.1906
#> s.e. 0.3395 0.2617 0.2346 0.1698 0.3816
#>
#> sigma^2 estimated as 3257345416: log likelihood=-418.67
#> AIC=849.34 AICc=852.45 BIC=858.5
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE
#> Training set 6516.114 51951.36 30469.81 -1.359721 23.11542 0.7101932
#> ACF1
#> Training set -0.03541161
OPSIONAL: Custom function compare_fitting() untuk membandingkan performa di data train agar tampilan lebih rapi
# opsional: custom function
compare_fitting <- function(x){
lapply(x, function(x) summary(x)["Training set", ]) %>%
lapply(., function(x) x %>% t() %>% as.data.frame) %>%
bind_rows() %>%
mutate(model = names(x)) %>%
select(model, everything())
}model_list <- list(
"ARIMA 011" = no_auto,
"ARIMA 111" = no_arima1,
"ARIMA 211" = no_arima2,
"ARIMA 411" = no_arima3)
compare_fitting(model_list)#> Series: no_train
#> ARIMA(0,1,1)
#>
#> Coefficients:
#> ma1
#> -0.8150
#> s.e. 0.0919
#>
#> sigma^2 estimated as 3317659997: log likelihood=-420.96
#> AIC=845.93 AICc=846.32 BIC=848.98
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE ACF1
#> Training set 10972.01 55929.24 31826.05 0.950838 23.63757 0.7418046 -0.04361855
#> Series: no_train
#> ARIMA(1,1,1)
#>
#> Coefficients:
#> ar1 ma1
#> 0.0267 -0.8216
#> s.e. 0.1943 0.1019
#>
#> sigma^2 estimated as 3420512259: log likelihood=-420.96
#> AIC=847.91 AICc=848.71 BIC=852.49
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE ACF1
#> Training set 11073.59 55922.5 31760.79 1.07617 23.49991 0.7402834 -0.06073625
#> Series: no_train
#> ARIMA(2,1,1)
#>
#> Coefficients:
#> ar1 ar2 ma1
#> -0.0445 -0.2210 -0.7471
#> s.e. 0.2000 0.1822 0.1416
#>
#> sigma^2 estimated as 3379874911: log likelihood=-420.28
#> AIC=848.55 AICc=849.93 BIC=854.66
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE
#> Training set 10008.72 54713.83 30326.12 0.6124584 22.41164 0.7068441
#> ACF1
#> Training set -0.06979143
#> Series: no_train
#> ARIMA(4,1,1)
#>
#> Coefficients:
#> ar1 ar2 ar3 ar4 ma1
#> -0.5841 -0.6503 -0.4078 -0.3948 -0.1906
#> s.e. 0.3395 0.2617 0.2346 0.1698 0.3816
#>
#> sigma^2 estimated as 3257345416: log likelihood=-418.67
#> AIC=849.34 AICc=852.45 BIC=858.5
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE
#> Training set 6516.114 51951.36 30469.81 -1.359721 23.11542 0.7101932
#> ACF1
#> Training set -0.03541161
Forecasting
Lakukan forecasting untuk semua kombinasi model
no_auto_f <- forecast(no_auto, h = 8) # panjang data test
no_arima1_f <- forecast(no_arima1, h = 8)
no_arima2_f <- forecast(no_arima2, h = 8)
no_arima3_f <- forecast(no_arima3, h = 8)Evaluate model
OPSIONAL: Custom function compare_forecast() untuk membandingkan performa di data test dengan tampilan lebih rapi
# opsional: custom function
compare_forecast <- function(x, test){
lapply(x, function(x) forecast::accuracy(x, test)["Test set", ]) %>%
lapply(., function(x) x %>% t() %>% as.data.frame) %>%
bind_rows() %>%
mutate(model = names(x)) %>%
select(model, everything())
}forecast_list <- list(
"ARIMA 011" = no_auto_f,
"ARIMA 111" = no_arima1_f,
"ARIMA 211" = no_arima2_f,
"ARIMA 411" = no_arima3_f)
compare_forecast(forecast_list, no_test)Visualisasikan hasil forecasting model
no_ts %>%
autoplot(series = "Actual") +
autolayer(no_auto_f$mean, series = "ARIMA(0,1,1)") +
autolayer(no_arima1_f$mean, series = "ARIMA(1,1,1)") +
autolayer(no_arima2_f$mean, series = "ARIMA(2,1,1)") +
autolayer(no_arima3_f$mean, series = "ARIMA(4,1,1)")Kesimpulan:
ts()adf.test(), harapannya p-value < alpha agar data stasionerauto.arima(obj_ts, seasonal=F)DAY 4
Seasonal ARIMA adalah metode ARIMA di mana data time series memiliki pola seasonal. Tahapan dalam melakukan pemodelan menggunakan SARIMA sama seperti saat membuat pemodelan ARIMA.
Keterangan:
Jika kita menggunakan model SARIMA, diperlukan langkah tambahan ketika melakukan differencing, yaitu differencing sesuai pola berulang/frequency untuk menghilangkan seasonality.
Data DAUTONSA.csv merupakan penjualan motor di suatu dealer dari tahun 1967-2019.
DATE: tanggal saat dilakukan pencatatan penjualanDAUTONSA: rata-rata penjualan bulanansales <- read.csv("data_input/DAUTONSA.csv")
glimpse(sales)#> Rows: 102
#> Columns: 2
#> $ DATE <chr> "2011-01-01", "2011-02-01", "2011-03-01", "2011-04-01", "2011~
#> $ DAUTONSA <dbl> 258.576, 338.615, 445.315, 405.418, 362.214, 348.950, 329.517~
Cek apakah data sudah memenuhi time series yang baik?
sales %>%
mutate(DATE = ymd(DATE)) %>% # ubah jadi tipe data date
arrange(DATE) %>% # mengurutkan berdasarkan tanggal
pad() %>% # menambal tanggal (padding)
anyNA()#> [1] FALSE
Artinya data kita sudah memenuhi karakteristik data time series yang baik
Data wrangling
DATE menjadi tipe data tanggalDATE untuk kebutuhan visualisasisales <- sales %>%
mutate(DATE = ymd(DATE),
MONTH = month(DATE, label = T))Visualisasi
sales %>%
ggplot(aes(DATE, DAUTONSA)) +
geom_point() +
geom_point(sales %>% filter(MONTH == "Jan"),
mapping = aes(DATE, DAUTONSA), col = "blue") +
geom_line()Membuat object time series
sales_ts <- ts(data = sales$DAUTONSA,
start = 2011,
frequency = 12) # data bulanan, pola tahunanDecompose
sales_ts %>%
decompose() %>%
autoplot()Cross-validation
# 1.5 tahun untuk test
sales_test <- tail(sales_ts, 1.5*12)
# 7 tahun untuk train, sisanya
sales_train <- head(sales_ts, -1.5*12)Cek stationarity data
Gunakan adf.test(): - H0 = Tidak Stationary - H1 = Stationary
sales_train %>%
adf.test()#>
#> Augmented Dickey-Fuller Test
#>
#> data: .
#> Dickey-Fuller = -2.9359, Lag order = 4, p-value = 0.192
#> alternative hypothesis: stationary
Kesimpulan: p-value > alpha (gagal tolak H0), artinya data sales_train masih tidak stasioner
ARIMA: diff(lag=1)
SARIMA: diff(lag=frequency)
Kasus frequency = 12:
sales_train_diff <- sales_train %>%
diff(lag=12) %>% # differencing untuk komponen seasonal
diff(lag=1) # differencing untuk trend
sales_train_diff %>%
adf.test()#>
#> Augmented Dickey-Fuller Test
#>
#> data: .
#> Dickey-Fuller = -4.7557, Lag order = 4, p-value = 0.01
#> alternative hypothesis: stationary
Kesimpulan:
Note: lakukan differencing seasonal terlebih dahulu, baru trend-nya agar komponen waktu tidak bergeser. Urutannya:
Baseline model
Gunakan auto.arima(), namun set parameter seasonal = T agar menjadi auto SARIMA
ARIMA: seasonal = F SARIMA: seasonal = T
sales_auto <- auto.arima(sales_train, seasonal = T)
sales_auto#> Series: sales_train
#> ARIMA(0,1,1)(0,1,0)[12]
#>
#> Coefficients:
#> ma1
#> -0.6468
#> s.e. 0.0900
#>
#> sigma^2 estimated as 892.7: log likelihood=-341.71
#> AIC=687.41 AICc=687.59 BIC=691.94
Auto SARIMA: ARIMA(0,1,1)(0,1,0)[12]
Visualisasi hasil fitting model
sales_train %>%
autoplot(series = "Actual") +
autolayer(sales_auto$fitted, series = "ARIMA(0,1,1)(0,1,0)[12]")Manual SARIMA
Gunakan custom function plot_cf() untuk menentukan order p, q, P, Q
# visualisasi PACF ACF
plot_cf(sales_train_diff, lag.max = 50)PACF
ACF
Kombinasi model ARIMA(p,d,q)(P,D,Q)[m]:
Fitting model manual
Parameter:
order = c(p,d,q)seasonal = c(P,D,Q)sales_sarima1 <- Arima(sales_train,
order = c(1,1,1),
seasonal = c(0,1,0))Gunakan custom function compare_fitting() untuk membandingkan performa model di data train
model_list <- list(
"ARIMA(0,1,1)(0,1,0)[12]" = sales_auto,
"ARIMA(1,1,1)(0,1,0)[12]" = sales_sarima1)
compare_fitting(model_list)#> Series: sales_train
#> ARIMA(0,1,1)(0,1,0)[12]
#>
#> Coefficients:
#> ma1
#> -0.6468
#> s.e. 0.0900
#>
#> sigma^2 estimated as 892.7: log likelihood=-341.71
#> AIC=687.41 AICc=687.59 BIC=691.94
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE ACF1
#> Training set -4.11363 27.27477 19.30697 -1.13262 4.466261 0.4645816 -0.02948427
#> Series: sales_train
#> ARIMA(1,1,1)(0,1,0)[12]
#>
#> Coefficients:
#> ar1 ma1
#> -0.0095 -0.6415
#> s.e. 0.1808 0.1361
#>
#> sigma^2 estimated as 905.6: log likelihood=-341.71
#> AIC=689.41 AICc=689.77 BIC=696.2
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE
#> Training set -4.092282 27.27413 19.30068 -1.128451 4.464672 0.4644302
#> ACF1
#> Training set -0.02586758
Forecasting
sales_auto_f <- forecast(sales_auto, h = 1.5*12)
sales_sarima1_f <- forecast(sales_sarima1, h = 1.5*12)Evaluasi model
Gunakan custom function compare_forecast() untuk membandingkan performa model di data test
forecast_list <- list(
"ARIMA(0,1,1)(0,1,0)[12]" = sales_auto_f,
"ARIMA(1,1,1)(0,1,0)[12]" = sales_sarima1_f)
compare_forecast(forecast_list, sales_test)Visualisasi hasil forecasting
sales_ts %>%
autoplot(series = "Actual") +
autolayer(sales_auto_f$mean, series = "ARIMA(0,1,1)(0,1,0)[12]") +
autolayer(sales_sarima1_f$mean, series = "ARIMA(1,1,1)(0,1,0)[12]")Optional: untuk visualisasi hasil forecast yang interaktif menggunakan fungsi test_forecast() dari package TSstudio
test_forecast(actual = sales_ts,
forecast.obj = sales_auto_f,
train = sales_train,
test = sales_test)Asumsi pada time series diujikan untuk mengukur apakah residual yang peroleh dari hasil modeling sudah cukup baik untuk menggambarkan dan menangkap informasi pada data. Mengapa menggunakan residual data? Karena kita dapat mendapatkan informasi dari data aktual maupun dari hasil prediksi menggunakan model.
Metode forecasting yang baik menghasilkan nilai residual berikut ini:
Untuk memastikan bahwa residual yang dihasilkan memiliki kriteria diatas, maka terdapat asumsi untuk mengujinya.
Note: Cek asumsi cukup dilakukan pada model akhir yang akan dipilih
Untuk cek ada/tidaknya autokorelasi pada hasil forecasting dapat menggunakan visualisasi plot residual ACF dengan fungsi Acf(model$residuals)
# $residuals selisih nilai actual dan forecast di data train
Acf(sales_sarima1$residuals)Berdasarkan plot ACF di atas, pada lag 8 dan 14 terdapat autokorelasi. Namun apakah secara keseluruhan lag, autokorelasi tersebut signifikan? Visualisasi bersifat subjektif, apabila ingin lebih pasti, lakukan uji Ljung-box menggunakan fungsi Box.test(model$residuals, type = "Ljung-Box") dengan hipotesis:
Harapan: p-value > alpha (0.05)
Box.test(sales_sarima1$residuals, type = "Ljung-Box")#>
#> Box-Ljung test
#>
#> data: sales_sarima1$residuals
#> X-squared = 0.058239, df = 1, p-value = 0.8093
Kesimpulan: p-value > alpha maka gagal tolak H0, sehingga tidak terdapat autokorelasi pada residual
Uji ini sama persis dengan uji normalitas pada course Regression Model.
hist(sales_sarima1$residuals, breaks = 15)Berdasarkan histogram di atas, distribusi dari residual masih rancu apakah normal atau tidak. Apabila ingin lebih pasti, lakukan uji shapiro test menggunakan fungsi shapiro.test(model$residuals) dengan hipotesis:
Harapan: p-value > alpha (0.05)
shapiro.test(sales_sarima1$residuals)#>
#> Shapiro-Wilk normality test
#>
#> data: sales_sarima1$residuals
#> W = 0.96604, p-value = 0.02575
Kesimpulan: p-value < alpha sehingga tolak H0, maka residual tidak menyebar normal (tidak memenuhi asumsi)
Keep in mind: model yang kita bangun masih bias
STLM = Seasonal Trend with Loess Model
Apabila dalam decompose biasa, dalam mendapatkan komponen trend dengan cara central moving average (CMA) di mana secara konsep setiap data yang ingin dirata-ratakan diberikan bobot yang sama sesuai ordo yang ditetapkan. Karena merata-ratakan data tengahnya hasilnya kita kehilangan data awal dan data akhirnya, sehingga ada beberapa informasi yang hilang.
Ada salah satu cara untuk mendapatkan decompose data namun tetap mempertahankan informasi dari seluruh data yang kita miliki yaitu dengan menggunakan STL (Seasonal Trend with Loess). STL secara konsep akan melakukan smoothing terhadap data tetangga setiap masing-masing observasi dengan memberikan bobot yang lebih berat terhadap data yang dekat dengan observed data.
Kekurangan dari STL hanya bisa melakukan decompose pada additive data, apabila terdapat multiplicative data dapat menggunakan transformasi log().
Untuk memodelkan hasil STL, kita juga dapat menerapkan metode exponential smoothing (ETS) dan ARIMA. Selain itu, STLM dapat digunakan sebagai alternative cara untuk menangkap seasonal yang belum bisa ditangkap oleh metode ETS dan ARIMA biasa.
Parameter stlm():
y : object time seriess.window : pola seasonal yang ingin ditangkap (frequency)method : metode forecast yang akan digunakan, tersedia ets dan arimasales_stlm <- stlm(y = sales_train,
s.window = 12,
method = "arima")
summary(sales_stlm$model)#> Series: x
#> ARIMA(2,1,1)
#>
#> Coefficients:
#> ar1 ar2 ma1
#> -1.3573 -0.5765 0.970
#> s.e. 0.0903 0.0897 0.046
#>
#> sigma^2 estimated as 562.9: log likelihood=-380.07
#> AIC=768.15 AICc=768.66 BIC=777.82
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE
#> Training set 0.4646435 23.15404 18.12287 -0.08624959 4.369887 0.4372299
#> ACF1
#> Training set -0.01882036
Berikut adalah hasil decompose dibalik STLM:
# decompose biasa
sales_train %>%
decompose() %>%
autoplot()# decompose stl
sales_stlm$stl %>%
autoplot()Forecasting
sales_stlm_forecast <- forecast(sales_stlm, h = 18)
forecast::accuracy(sales_stlm_forecast$mean, sales_test)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -37.34721 46.99031 38.06995 -11.28805 11.52644 0.4927543 0.9266351
Visualisasi
sales_ts %>%
autoplot(series = "Actual") +
autolayer(sales_stlm_forecast$fitted, series = "Train") +
autolayer(sales_stlm_forecast$mean, series = "Test")ARIMA(p,d,q)(P,D,Q)[m]
auto.arima(obj_ts, seasonal=T)